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Abstract. We study the onset of the reheating epoch at the end of axion-driven inflation 
where the axion is coupled to an Abelian, 17(1), gauge field via a Chern-Simons interaction 
term. We focus primarily on m 2 (j) 2 inflation and explore the possibility that preheating 
can occur for a range of coupling values consistent with recent observations and bounds on 
the overproduction of primordial black holes. We find that for a wide range of parameters 
preheating is efficient. In certain cases the inflaton transfers all of its energy to the gauge fields 
within a few oscillations. In most cases, we find that the gauge fields on sub-horizon scales 
end preheating in an unpolarized state due to the existence of strong rescattering between the 
inflaton and gauge-field modes. We also present a preliminary study of an axion monodromy 
model coupled to 17(1) gauge fields, seeing a similarly efficient preheating behavior as well as 
indications that the coupling strength has an effect on the creation of oscillons. 
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1 Introduction 

Inflationary model building has entered a particularly exciting phase with the demonstra¬ 
tion by the BICEP2 experiment [1] of the sensitivity to B-mode polarization of the cosmic 
microwave background (CMB) at a level where interesting constraints can be, or soon will 
be, placed on the inflationary model space. The recent Planck results on dust emission [2] 
combined with the joint Planck and BICEP2/KECK Array analysis [3] mean that dust- 
foregrounds need to be accurately characterized in order to determine if any of the B-mode 
polarization observed by the BICEP2 experiment is due to primordial gravitational waves. 
If confirmed, the observation of these gravitational waves present a spectacular confirmation 
of one of the early observational predictions for slow-roll inflation [4], However, they also 
present challenges for inflationary model building. The large primordial gravitational wave 
amplitudes required to explain the BICEP2 signal generically require that the scalar-inflaton 
field rolls over a distance in field space that is large compared to the four-dimensional Planck 
scale [5]. This makes the theory extremely sensitive to unknown physics at short wavelengths, 
e.g. in the ultra-violet (UV), at or near the Planck scale, and leads to a loss of predictive 
power, if not a loss of the inflationary mechanism itself. A possible way around this UV 
sensitivity problem is to find a symmetry powerful enough to forbid interactions between the 
sector driving the inflationary expansion and other unknown physics. 

A promising candidate for the inflaton field is a pseudo-scalar or axiom These fields 
enjoy shift-symmetries that protect their role as inflatons from being spoiled by coupling to 
unknown UV physics. Shift-symmetries require that the theory is invariant under a constant 
shift of the field value, and severely restrict the form of possible interactions with other 
fields. One of the earliest proposed axion inflation models was Natural inflation [6, 7]. In 
this scenario, a cosine potential for the axion is generated by the condensation of a non- 
Abelian gauge group. Slow-roll inflation is achieved by the hierarchy between the height and 
width of the potential, the form of which is protected by the shift-symmetry. In order to 
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generate density fluctuations with a spectrum that reproduces the observed fluctuations in 
the CMB, the axion is required to have a periodicity, or associated mass scale, larger than 
the Planck scale. This super-Planckian periodicity makes it extremely difficult to embed 
Natural inflation in a fundamental theory such as string theory [8]. Model builders sought to 
circumvent this obstruction in a model dubbed IV-flation [9-11] by using a large number of 
axions, each with small periodicities, to implement assisted inflation [12]. While each of the 
N axions only rolls a small distance in field space, their collective motion is responsible for 
inflation and effectively traverses a large distance. Unfortunately the N axions in IV-flation 
have the effect of renormalizing the Planck scale to a lower value, and the theory ultimately 
suffers from similar pathologies as the original Natural inflation model it was supposed to 
fix. Other formulations making use of misaligned axions have also been proposed [13-15]. 
More recently, the observation that a single axion undergoing monodromy can have a small 
periodicity while ultimately traversing a large field range [16, 17] has seen resurgent interest in 
axion inflation [18-23]. For a recent review of axion inflation see Ref. [24] and its realizations 
in string theory see Ref. [25]. 

At the end of the inflationary phase, the Universe must undergo a phase transition, from 
its super-cooled state to a state filled with radiation and ultra-relativistic matter, to begin 
the hot big bang. The physics of this phase transition is thought to be highly non-linear, 
and its details are unknown (see e.g. Ref. [26] for a recent review). The shift-symmetry in 
axion-driven inflation that is so effective at protecting the form of the inflationary sector 
from unknown UV physics also severely constrains the form of its couplings to the visible 
sector. These couplings to the standard model of particle physics, either directly or indirectly 
via intermediaries, are required in order to transfer the inflaton energy into radiation and 
ultra-relativistic matter. The shift symmetry dictates that couplings to matter fields must 
be derivative interactions, therefore a class of allowed interactions are those in which the 
axion is coupled to a gauge field through a dimension-five Chern-Simons term or Pontryagin 
density. A coupling of this form is allowed by the symmetries and, from the viewpoint of an 
effective field theory, must be present. Furthermore, this coupling provides a perturbative 
decay channel for the axion into gauge bosons which guarantees that reheating eventually 
completes through perturbative decays alone, and thus provides a viable pathway through 
which the Universe can transition from inflation into the hot big bang. 

The effects of the coupling of axions to gauge fields during inflation is, by now, a well 
studied field. It has been known for a long time that axion-gauge field couplings lead to parity- 
violating gauge-fields that are amplified during slow-roll inflation [27, 28]. The behavior 
of these gauge fields during inflation and their influence on the inflationary dynamics has 
also been extensively studied [29-34], Further, the authors of Refs. [35, 36] studied metric 
fluctuations generated by a rolling auxiliary pseudo-scalar during inflation. Some work has 
also been done on the reheating of axion-driven inflation. The authors of Refs. [37] considered 
stringy reheating of a monodromy scenario while Ref. [21] studied perturbative decay of 
an axion inflaton to photons/gauge bosons. Furthermore, previous studies have considered 
perturbative analyses of the Mathieu equation for U( 1) gauge fields coupled to an oscillating 
scalar or collection of scalars [38, 39]. However, these studies were focussed primarily on 
ranges of parameters considered natural for scenarios of Natural inflation and A-flation and 
concluded that highly non-linear effects and parametric resonance were unimportant in these 
models. 

Recently, larger couplings between the axion and gauge sectors have been considered 
[30, 32], which can lead to observable effects during inflation due to the rescattering of the 
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gauge fields off the axion condensate [40]. Effects such as non-Gaussianity of the density 
fluctuations, chiral gravitational waves, and the production of primordial black holes [31, 32, 
40-42] place upper-bounds on the strength of the couplings, the most stringent being due to 
restrictions on the production of primordial black holes at the end of inflation. 

Our current work addresses a specific problem; given that high-scale inflation requires 
highly constrained couplings, can the Universe undergo a preheating phase (either tachyonic 
or resonant) following axion inflation? Numerical investigations of reheating in canonical [43— 
55] or non-canonical [56] scalar field scenarios is becoming routine. Furthermore, couplings 
of scalar fields to Abelian [57-59] or non-Abelian gauge fields in cosmological settings have 
been studied in recent years [60-75]. The question we address here is new. We employ 
lattice simulations, using the same numerical technique of Ref. [59], to study the possibility 
of tachyonic or parametric amplification of gauge fields after inflation due to the presence of 
a Chern-Simons term. As in most studies of preheating, we remain agnostic as to the specific 
nature of the coupled gauge field and do not necessarily identify it as a 17(1) gauge field from 
the standard model. 

Our results can be summarized as follows. We find that for reasonable ranges of the 
axion-gauge field coupling, non-linear effects can be very important at the end of inflation. 
In particular, at the middle to upper range of the couplings allowed by black hole abundance, 
we find that reheating is essentially instantaneous, proceeding via a phase of tachyonic reso¬ 
nance [76] and completing within a single oscillation of the axion. Despite the asymmetry in 
the equations of motion for the two polarizations of the gauge fields, on sub-horizon scales, 
rescattering of the gauge bosons off the axion condensate is efficient at generating the second 
polarization. On scales larger than the horizon at the end of inflation, an asymmetry between 
the gauge field polarizations remains. The Universe that results in these cases is radiation 
dominated and is characterized by a very high reheating temperature. As the coupling is 
decreased, this tachyonic resonance is weakened and the axion oscillates multiple times be¬ 
fore reheating completes. During these multiple oscillations equal levels of both polarizations 
of the gauge field are excited. Decreasing the coupling further yields a brief window where 
parametric resonance effects become important before preheating abruptly shuts off and non¬ 
linear effects cease to be important. At these lower couplings, non-linear effects are negligible 
and the Universe reheats via perturbative decay of the axion into gauge bosons. 

We also investigate how these preheating effects might depend on the shape of the 
inflationary potential. As a second test, we subject the axion to a monodromy-type potential 
[16], and show that the range of couplings for which efficient preheating can occur is slightly 
different, although the same order-of-magnitude. We conclude by presenting an intriguing set 
of data that suggest gauge fields might play a role in the creation and evaporation of oscillons 
in this scenario. 

We work in natural units where h = c = 1, however, we retain the Planck mass, m p \ = 
1.22 x 10 19 GeV. 


2 Background and conventions 

We begin with the usual action for axion-driven inflation 


Sinf = J d A Xy/^g 


m 


. i 
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where 0 is a pseudo-scalar (axion) and V((j>) is a potential that supports slow-roll inflation. 
For definiteness, we consider the potentials for the simplest type of chaotic inflation [77], 

V{ct>) = *m 2 0 2 , (2.2) 

and the simplest type of axion monodromy inflation [16] 

v (</>) =/i 3 ( V4> 2 + 4>l- 0c) , (2.3) 


which is well described by a linear function of (j) for large field values. The amplitude of the 
scalar spectrum fixes the parameters m and /i to be [78-80] 

m « 1.06 x 10 -6 m p i, (2.4) 

and 

H ~ 1.20 x 10 _4 m p i. (2.5) 

The parameter (f) c in the monodromy potential, Eq. 2.3, has a negligible effect on the spectrum 
of curvature fluctuations and gravitational waves on the scales that are observable in the CMB 
(provided, of course, that 0 C is much smaller than the held value where the CMB fluctuations 
are generated) and, consequently, is unconstrained by data. However, (j) c becomes important 
near the end of inflation and has a small effect on the value of (j) when inflation ends. Further, 
for small held values, (j) < (j) c , the potential can be expanded as 
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and the resulting dynamics of 4> in this region depend strongly on the value of (j) c [81]. 
In addition to the axion, we consider a 17(1) gauge held coupled to the axion 


^gauge — j 7 X\J Jj 


^iF^ 


(2.7) 


where a is a dimensionless coupling constant of order unity and / is a mass scale associated 
with the pseudo-scalar (axion). The held strength and its dual are given by the standard 
expressions, F /IV = d^A v — and F^ F a ^/2, where is the completely 

antisymmetric tensor, and our convention is 


0123 _ 1 

" y/=g‘ 


( 2 . 8 ) 


Greek letters here and throughout denote four dimensional spacetime indices and Roman 
letters from the middle of the alphabet are used to denote spatial indices. Repeated lower 
spatial indices are summed using the Kronecker delta. We work with the Friedmann-Lemaitre- 
Robertson-Walker (FLRW) metric in conformal time with mostly-plus conventions 


ds 2 = —a 2 (dr 2 — dx 2 ). 


(2.9) 


The equation of motion for the pseudo-scalar held is the Klein-Gordon equation sourced by 
the Chern-Simons density of the gauge held 


(■ d 2 + 2 Hd T 


didi)(j) + a 2 


dV 

d<j> 


a 

47 


a 2 F^F^, 


( 2 . 10 ) 
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where % = a'/a and where here and below ' = d T = d/dr. The equations of motion for the 
gauge field are 


_ r\ _ ^ 

d p + -d p {^c^Fn = o. 


The cr = 0 equation is the Gauss’ law constraint 


a 


djdjA 0 - d T diAi + -e^d^diA^ = 0 , 


( 2 . 11 ) 
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while the a = i equations are the field equations for the spatial components of the gauge field 

Q> Oi 

8t (drAi 5jj4o) T 8m(8m-Ai 9jA m ) + — (imk^r^dmAk ~J ^imk^m4 > i.^T-^k ^fcV^o) —0. 

(2.13) 

Finally, assuming the metric is unperturbed, the scale factor satisfies Einstein’s equations 


3 777,: 


87T 


pl % 2 = a 2 p, 


and 




! P + P 


8vr v ' 2 

The pressure, p. and energy density, p, are found from the stress-energy tensor 


Tpv =Tr [Fp a F u p\ g ap - 9 -fF p „F^ - g^ 


which can be explicitly written as 
1 <\>' 2 , 1 (< 9 i ^) 2 


T>g prT d p (l)d a 4> + v{(p) 


and 


P ~2^ + 2^t~ + V ^ + 2^ {d ° Ai ~ diAo) 2 + 4 ^ 9iAj ~ djAi)2 ' 


+ v w + - **> 2 + 1 ir^ - ^ )2 - 


(2.14) 

(2.15) 


+ dpcjidvd, (2.16) 


(2.17) 


(2.18) 


Note that the axion-gauge field coupling does not contribute directly to the stress-energy 
tensor. 


3 Gauge-field production during inflation 

During inflation, the coupling of the gauge field to the axion results in exponential production 
of one polarization of the gauge field over the other [27-29]. To see this, we first fix the gauge 
by choosing Coulomb, or transverse, gauge d{Ai = 0. The Gauss’ law constraint, Eq. 2.12, 
then implies that Aq = 0 at linear order in fluctuations. 1 

1 Note that at this order (linear) in fluctuations, Coulomb gauge and temporal gauge (Ao = 0) are equiva¬ 
lent. This can be seen trivially from Eqn. (2.12), the Gauss’ law constraint. In the linear regime, this equation 
reads djdjAo — d T djAj = 0. In Coulomb gauge, assuming k ^ 0, this constraint reads Ao = 0. In temporal 
gauge, Gauss’s law reads d T djAj = 0, which for k ^ 0, implies djAj = 0. 
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At linear order in fluctuations, with this choice of gauge, the equation of motion for the 
gauge field becomes 


a 


d T Ai O rn () m Aj ^ timk9 T 0^rn A/,- — 0 


(3.1) 


To study the fluctuations of the gauge field, it proves most convenient to work in Fourier 
space, where our convention is 


A( X ) = 


d 3 k 

(2ir) 2 


A k e 


i k x 


The Fourier components are then expanded in a basis of helicity states 

A k = ^ A^s A (k), 


(3.2) 


(3.3) 


A=± 


where the polarization vectors e A (k) satisfy the orthogonality and normalization relations 

hefi k) =0, e ljk kj£^( k) = ^fikef( k), 

£ t ( k )* = £ t (- k )> 4( k ) e * V (- k ) = (3-4) 

In this situation, conformal time is defined to be a negative, increasing quantity during 
inflation 


dr = 


dt 


' dt 


t = 


din a 
aH 


1 

aH' 


(3.5) 


where the last approximation is exact in the de Sitter limit, eh —> 0, where the slow-roll 
parameter, e# is defined as en = —H/H 2 . Here and throughout, an overdot is used to 
denote a derivative with respect to cosmic time, t. 

We can now quantize the modes by introducing the creation and annihilation operators, 
a A (k) and at (k) satisfying the canonical commutation relations 


a A(k),a^,(k') = (27r) 3 5 AA /d 3 (k - k'), 

which allows us to expand the mode-functions as 


Ai(r, x) = 


A=± • 


d 3 k 

(27r) 3 


e lk ' x e A (k) A A (k, r)a A (k) + A A ’*(k, r)a^(-k) 


A,* / 


(3.6) 


(3.7) 


With our conventions, the gauge field equation of motion Eq. 2.13 becomes a separate equation 
for each polarization, depending only on the magnitude of the momenta k = |k 


^2 , , 2 ■ a </> k 


dr + k^j— A ^ = 0 , 


(3.8) 


where we have also made use of the de Sitter approximation for the scale factor during inflation 
in the last term. During inflation, en = 47Ti^ 2 /(i7 2 mp|) ~ const, and, after changing variable 
to u = 2 ikr, the equation of motion is transformed to the Whittaker equation 

f d\ 1 | A | 1/4 — fi 2 


\dz 2 


z 


W X!tl (z) = 0 . 


(3.9) 










In our case, we have 


4 ^ l u ) 


and thus fi = 1/2 and A = q=i£ where we have dehned 

la 4> . rripi a 

‘=2jH =s ' snW 7m7 



(3.10) 


(3.11) 


The general solution of the Whittaker equation can be written in terms of the Whittaker 
W-function 


Af(kr) = C 1 W^i(2ikr) + C 2 W ±i( , i (-2ikr). (3.12) 

The constants of integration, C\ and C 2 , are set by canonical quantization which amounts to 
normalizing the modes according to the Wronskian condition 

W[A±(kT),(A±(kT))*]=i, (3.13) 


and demanding that the modefunction approaches the Minkowski vacuum in the limit, k\r\ —> 

oo 


lim A ± (k, r) = —jL= 

fc|r|—>oo '/2k 


D —ikr^fi^ ln(—2 kr) 


The properly normalized solutions are 


e 2 ^ 


A ( k ’ T ) = ~^ w Ti^ 2ikT ^- 

In the limit k\r\ —>• 0, the asymptotic form of the mode function is 

, 1 e ±f* 

lim A ± {k,r) = ,_ ———■——. 

fc|r|->-o v^r(i±*o 


(3.14) 


(3.15) 


(3.16) 


Compared to the conformally invariant radiation solution, the circularly polarized modes get 
amplified by a factor 


A ± 

^4±,rad 


~ e2 


l€l±f£ 


(3.17) 


where we have used the Stirling formula 

|T(1 =F *01 - (2vr|^|) 1/2 e _7r| ^ l/2 , (3.18) 


and we have assumed that |£| > 1. Note that this means that when £ > 0 (£ < 0), the mode 
At {AT) gets amplihed by a factor ~ e 7r ^l while the other mode is unchanged. For this work, 
we focus on large-held inhationary models, and assume that (j) < 0 so that the mode A~ is 
amplihed during inflation. 

The exponentially enhanced gauge helds have important effects during inflation due to 
their re-scattering off the inhaton condensate and their interactions with the metric. The for¬ 
mer leads to the production of fluctuations of the inhaton which are statistically non-Gaussian, 
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while the latter leads to the production of gravitational radiation [32], The Gaussianity of 
the observed density fluctuations by Planck [82] then implies that the quantity £cmb 2.22 
[24], where £cmb is the quantity in Eq. 3.11 evaluated during the time when the modes that 
form the CMB leave the horizon. 

During inflation, en and thus the ratio <f>/(Hm p i) increases. This means that shorter- 
wavelength modes that leave the horizon later in inflation are amplified more than their 
longer-wavelength counterparts that leave the horizon earlier. The largest effects occur when 
€h is near unity near the end of inflation. In the limit that £ 1 (which for models that 

satisfy £cmb <2.22 only possibly occurs near the end of inflation) the energy density in the 
gauge fields becomes important and the gauge-held fluctuations begin to backreact on the 
homogeneous background equations of motion. In this limit, in the Hartree approximation, 
the Friedmann (Eq. 2.14) and Klein-Gordon equations (Eq. 2.10) become 

+ ^ + (3.19) 

o7r z z 

-n 2 )=-[^- + ‘ 1 -a 2 {E 2 + B 2 )^, (3.20) 

and 

cj)" + 2n^' + a 2 V' = ja 2 (E- B), (3.21) 

where the electric and magnetic fields 2 associated with the 17(1) gauge held are E{ = a _2 7l( 
and Bi = a~ 2 eijkdjAk. In this limit, up to an irrelevant constant phase, the gauge held mode 
that is amplified is approximated near horizon crossing by [30] 

A k ( T > = 7^7 exp (vr|£| - 2 v / 2|C|fc|r|) , (3.22) 

while the other mode is unaffected and is negligible. The expectation values of the quantum 
helds are well approximated by [30] 3 

\{E 2 + B 2 ) ~1.4 • 10- 4 ^e 2 ^l, (E- B) ~ 2.4-10” 4 ^e 2 ^l. (3.23) 

Toward the end of inflation, for large values of m p \ot/ f, the back reaction of the gauge helds 
on the rolling axion becomes important and inflation is prolonged [32], During this phase, the 
primordial density fluctuation spectrum is expected to be dominated by rescattering and large, 
non-Gaussian density fluctuations are predicted. These large amplitude density fluctuations 
can produce primordial black-holes [83-85] and ensuring that they are not overproduced 
requires that £cmb 1-5 — 1.7 for m 2 q i 2 [41, 42] and £cmb 1-8 for monodromy [42], which 
is tighter than the current bounds from the Gaussianity of the CMB fluctuations. These 
limits are model dependent, and are somewhat sensitive to the form of the potential. In this 

2 We refer to these fields as electric and magnetic, however, they need not be the electro-magnetic fields of 
the standard model. 

3 Note that the axion velocity assumed here is opposite to that assumed by [30]. This means relative to 

this work, the other gauge mode is amplified and consequently the sign of (E ■ B) is opposite. 
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work we consider only values of the coupling a/f such that £cmb < 1-5 — 1.7. This bound 
translates to roughly a/f < llOm^ — 125 m^ 1 for the m 2 (f 2 potential and a/f < 180 
for the simple monodromy potential of Eq. 2.3. In practice, we do not approach this threshold 
in our simulations, keeping the couplings used for our simulations lower by around a factor 
of two. 

The bounds from primordial black-hole abundances rely on the above approximations of 
Eqs. 3.19 - 3.21 together with Eq. 3.23 being a good description of the system in the region 
of strong back reaction [41]. However, these approximations do not yield a self consistent 
set of equations because the resulting stress-energy tensor is not covariantly conserved. This 
means that in the region where back reaction becomes significant, the above equations are 
inaccurate. For small couplings, we expect that the above approximations are accurate enough 
to capture the onset of back reaction. In what follows we initialize our lattice simulations 
using the values of the fields found from the numerical evolutions of the Eqs. 3.19 - 3.21. 
We make use of the approximations of Eq. 3.23 for the expectation values of the energy 
density in gauge field fluctuations and the Pontryagin density respectively. The use of these 
approximations implies that the initialization of our simulations is less and less accurate for 
larger values of the coupling between the gauge field and the axiom To minimize this error, we 
initialize our simulations two e-foldings before inflation ends. The modes that are important 
in the reheating era are well within the horizon at this time, and well described by the linear 
approximation at the start of our simulations (See Table 1 and Fig. 6). 


4 Perturbative reheating 


Before we move to the non-linear regime, we briefly revisit the perturbative reheating case. 
The coupling of the axion to gauge fields provides a natural decay channel for the axion 
to produce gauge bosons. Even in the absence of non-linear effects due to the homogenous 
motion of the inflaton condensate, the Universe will reheat to gauge bosons via perturbative 
decays of the inflaton into a pair of gauge bosons. When the Hubble rate becomes comparable 
to the perturbative decay rate, the decay of the inflaton into gauge bosons proceeds rapidly. 
The contribution of these newly created particles to the energy density quickly dominates 
and reheats the universe. This represents a lower bound on the temperature of the Universe 
at reheating. 

The (j)FF coupling allows the axion to decay to two gauge bosons with a rate that is 
well known (see e.g. [86]) 

2 S 

am, 

r 4-+aa - 

where m</, is the mass of the axion about its vacuum. This rate sets a lower bound on 
the reheating temperature which is found by comparing the decay rate to the Hubble rate; 
perturbative reheating finishes when T/3i7 ~ 1, which results in a reheating temperature 


Treh ~ 



(4.2) 


where g * is the number of relativistic degrees of freedom. For the simplest version of chaotic 
inflation, we can write 


T reh ~ 1.31 x 10 9 






\ 3/2 


1.06 X 10 


-6 


ro-pl. 


«/(/ro p i)\ 


10 


GeV. 


(4.3) 
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For small values of the coupling, this is the dominant effect. Reheating in this case proceeds 
by perturbative decay of the axion into gauge bosons. 

5 Instabilities and resonance 

We are interested in what happens immediately following the inflationary epoch. We assume 
that the pseudo-scalar begins to oscillate about the minimum of its potential, which we 
initially take to be a quadratic function, Eq. 2.2. While this form is exact for the simplest 
model of inflation, there are anharmonic corrections to the potential for other important axion 
inflation scenarios, such as monodromy inflation. 

It is instructive to first consider the behavior of the system by temporarily neglecting 
the expansion of space and considering only the linear theory to gain intuition for the system 
and to identify features present in the fully nonlinear treatment. In this regime the axion 
satisfies the equation 

cj) + m 2 4> = 0, (5.1) 

where overdots denote derivatives with respect to cosmic time t. This equation has the simple 
solution 


cj>(t) = 4> o cos (mt), 


(5.2) 


where </>o is (approximately) the field value at which the slow-roll conditions are violated 
and inflation ends. To estimate the effect of parametric resonance (while still neglecting the 
expansion of space) we write the equation of motion for the mode amplitudes, Eq. 3.8, as 


Ai+k\k=f 


a 


*£ = o- 


(5.3) 


With the solution for <f>[t) from Eq. (5.2), and after redefining time z = mt/ 2, this equation 
can be recast as 


' d 2 4 _fc 
dz 2 m 


( k a 

— T T 0ocos(2z) 

\m f 


At = 0 . 


(5.4) 


The fact that each helicity obeys a different equation is irrelevant here because the sign of 
the term proportional to a/f in Eq. 5.4 can be reversed by a constant shift of its argument 
z — > z + 7t/ 2. In this approximation both modes are expected to grow equally after each 
inflaton oscillation. We show how this symmetry between the two helicities is broken once 
the expansion of the Universe is taken into account. 

We can compare Eq. 5.4 to the normal Mathieu equation, 

d 2 u , . ... 

^2 + [A k + 2 q cos(2z)] u = 0, 

to see that 4 

Ak = 4:(—'\ , <? = t2—^00- (5.6) 

\ rn J m j 

4 This definition of the two parameters of the Mathieu equation as A *, and q is typically used in models 
where the inflaton decays to scalars. In these models q does not depend on the wavenumber. In the case of 
gauge fields q does depend on k, but we refrain from using a subscript to be consistent with prior literature. 


(5.5) 
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From Eq. 5.5 we see two important thresholds that define the behavior of the solutions. 
These thresholds are the tachyonic resonance threshold, which is set by A & < 2 q, and the 
broad-to-narrow resonance threshold, which is determined by the size of q. Broad resonance 
occurs for q 1 while narrow resonance occurs for q < 1. Fig. 1 shows the Mathieu 
instability diagram for our process, along with three { q , A *,} curves for fixed coupling, a/ f(f> o, 
and varying wavenumber. The A = 2 q line is the diagonal for our choice of axes range. It 
can be clearly seen that in the regime where A & < 2 q, the instability bands are much broader 
and the imaginary parts of the Mathieu exponents are larger. It is also interesting to note 
that this system can be cast in terms of only two (dimensionless) combinations: the ratio of 
wavenumber to mass scale, k/m, and the product of the coupling strength and initial axion 
amplitude, («//)/>o- The curves in Fig. 1 are defined by a fixed coupling and initial held 
amplitude (a/f)cj) o, and are parameterized by wavenumber. 



Figure 1 . A Mathieu instability diagram for the process in question. The left panel is a contour 
plot of the magnitude of the imaginary part of the Mathieu characteristic exponent, Im(/r), which 
paramaterizes the rate of growth of the instability-non-zero values of this quantity indicate modes that 
grow exponentially. The three curves on the plot show the slice of parameter space along which the 
momentum modes in our system lie for three different couplings, a/f = 35 m“/ where <f> o ~ 0.20 m p i 
(black), a/f = 45 777 “/ where 4>o « 0.22 m p i (red), and a/f = 55771 “/ where fo ~ 0.24 to p i (blue). 
The right panel shows the size of the characteristic exponent as a function of the momentum of the 
mode k/m for the same three couplings. These curves identify the modes that are excited during 
the first oscillation of the held; as the amplitude of the oscillation decreases, these curves move to 
smaller q and become steeper. The range on the right panel is chosen to correspond to the values of 
the comoving wavenumbers that can be probed by our simulations. 


When the combination A + 2 q cos(2z) < 0, which occurs for A < 2 q, there is a 
tachyonic instability in the equation of motion [76]. This condition defines a set of modes 


k 

m 


< 7 <h ’ 


(5.7) 


whose mass-squared is negative (u" oc u from Eq. 5.5). There is always a set of such modes, 
as long as the homogeneous mode of the axion is oscillating. 

The details of this tachyonic regime are seen directly from the equations of motion for 
the two polarizations, Eq. 5.3. When the axion velocity changes sign, the combination 


cl ■ 

k^j(t>< 0 , 


(5.8) 
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Figure 2. The left panel shows the instability diagram for our process, reparameterized in terms of 
k/m and <poa/f. The blue, red and black lines correspond to the same couplings as in Fig. 1. The 
yellow curve divides the regions of narrow and broad resonance defined by q = 1 and the green line 
corresponds to (a/f)(f >o = 1, which is the value used in [38] . The right panel shows the instability 
chart in the region of small coupling (a/f)(f> o and small k/m. Note that the color scale is not the 
same in both plots. In the left plot, the largest value of the imaginary part of the Mathieu exponent 
is 3(/i) = 4.2 while in the right plot it is = 0.45. 


is negative for only one of the two polarizations, depending on the sign of <f>. When the axion 
velocity is positive (negative), the A + (A - ) mode is tachyonically amplified. This regime is 
obviously most important for large couplings and when the amplitude of the field oscillation 
is large. In these large coupling cases, preheating is extremely efficient and can complete 
after only one or two oscillations of the homogeneous inflation condensate. These tachyonic 
instabilities disappear when the homogeneous mode of the inflation breaks down, which occurs 
due to back scattering of the gauge modes, or self resonance of the axion itself (when the axion 
potential includes anharmonic terms). This should happen very quickly for large couplings 
- in some cases, only one polarization is amplified by this tachyonic resonance. However, 
as we show, in these cases rescattering or backscattering effects are extremely efficient and 
generate equal amounts of the non-tachyonic. polarization. As can be seen from Fig. 1, the 
wavenumbers corresponding to the tachyonic regime (Eq. 5.8) lead to much larger Mathieu 
exponents, so they dominate the behavior of the system, hence we concentrate on them. 

Having discussed the behavior of Eq. 5.5 for A < 2 q, we move to the second threshold 
which is defined by the size of q. If the {q,Ak} curve (with constant a/fcj) o) intersects the 
Mathieu bands for values of q 1, we are in the broad resonance regime [76, 87] characterized 
by the curve intersecting large instability bands. In the opposite regime, q < 1, the modes 
that are amplified have frequencies comparable to those of the oscillating inflaton. This is 
called the narrow resonance regime, since the instability bands of the Mathieu chart have a 
narrow width. The growth rate of gauge field modes in this regime can be fully analyzed using 
the methods of parametric resonance based on Floquet theory, as described for example in [87] 
and applied to narrow resonance for gauge fields in [38] . Since narrow parametric resonance is 
only present for small values of the coupling where the Universe does not completely preheat, 
we do not consider it for the remainder of this work. 

Previous studies of reheating after axion inflation inflation through gauge-field produc¬ 
tion, such as [38], focused on the region where q < 1, and seem to have missed the efficient, 
tachyonic preheating phenomenon on which we focus here. As can be seen in Fig. 2, both the 


- 12 















range of wavenumbers, as well as the size of the Floquet exponent are much larger for the 
range of couplings we study. 

These tachyonic instabilities exist for long wavelength modes of the gauge field, but the 
assumption of a homogeneous inflaton, 4>(t), breaks down quickly once back-reaction occurs. 
To see whether rapid production of gauge field quanta occurs, and/or whether this final state 
is polarized, quickly becomes a numerical question. However, there is some progress to be 
made using semi-analytic methods by reintroducing the expansion of the Universe into the 
equation of motion for the gauge fields, as explained in the next section. 

5.1 Semi-analytic treatment 

A detailed description of tachyonic resonance in the static-universe approximation can be 
found in [76], where it is shown that the agreement between these analytic results and nu¬ 
merical simulations (neglecting re-scattering and back reaction) is excellent (for g> 1 and 
Ak <2 q — 2y fq). In our current study, we go beyond this approximation. The richest phe¬ 
nomenology comes from the period between the end of inflation and the first few oscillations of 
the inflaton. During this epoch, the scale-factor is evolving nontrivially in time and cannot be 
approximated by an exponential, as it could during slow-roll inflation, or a power-law, as it is 
in a radiation-dominated universe after reheating. In order to proceed with a semi-analytical 
treatment, as an approximation, we model the evolution of the Universe during this time by 
the evolution of the classical axion in a background FRW spacetime, neglecting the effect 
of back reaction of the gauge fields on the expansion. In Fig. 3 we plot the evolution of the 
inflaton field, its velocity, and the Hubble parameter in units of the inflaton mass m for simple 
chaotic inflation. Without loss of generality, consistent with Section 3, we chose the value of 
the axion to be positive during inflation, (f> > 0, and its time derivative to be negative, (j> < 0 
(note that here VU > 0). This then determines the signs of these two quantities immediately 
at the end of inflation. Note that during the first oscillation the amplitude falls by about a 
factor of two, and thus approximating the inflaton background by a sinusoidal function with 
a constant amplitude is insufficient. 



0 2 4 6 8 10 


t (m 1 ) 

Figure 3. The evolution of the axion field 4>/m p i (blue solid), the velocity 4>/m p \ (red dashed) and 
the Hubble parameter H = a/a (black dotted) after the end of inflation, t = 0. Time and related 
parameters are measured in units of the axion mass to. 

We need to extend the method of [76] in order to treat the case of a time-dependent 
but non-harmonic effective frequency. We are most interested in the cases where the growth 
of the gauge field modes is due to tachyonic effects near the end of inflation. In these cases, 
preheating is extremely efficient, and we focus our attention on the growth of the gauge field 
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modes during the first oscillation of the axiom Due to the strong growth of the fluctuations, 
analytical methods quickly fail and we enter the regime of strong back-reaction, dominated 
by non-linear physics. In this regime we have to rely on full numerical simulations, which 
we present in the following section. Extending the (linear) analysis beyond the first axion 
oscillation is straightforward, and can be achieved by combining the procedure we describe 
here with the method of [76]. 

We start from Eq. 3.8 written in cosmic time 


At + HAt + 


and redefine the gauge held, = a l / 2 A± 


k\ a k ■ 

- Tv ~<t> 
a fa 


A ± = 0, 


Xk + 


k 2 a 2 a \ ak ■ 

W + 779 “ 7T7 ) ^ 7 :^ 


Xfc =o, 


a* 4 a 2 2aJ ' fa 
which leads to the equation of motion for this rescaled gauge held 

•■± , r,.±/ni 2 „± 


S+KW] Xfc=°- 


(5.9) 


(5.10) 


(5.11) 


We use the Wentzel-Kramers-Brillouin (WKB) approximation to describe the solution 
to Eq. 5.11 in the regions where the frequency 



f k 2 a 2 a \ a k ■ 
U 2 + 4a 2 2a) ^ f a ’ 


(5.12) 


is varying slowly, that is, |w/u; 2 | <C 1. Following the general procedure of the WKB ap¬ 
proximation, we distinguish three regimes based on the behavior of the effective frequency. 
We track a mode as it goes from having a positive frequency-squared, [w^(t)] 2 , (regime I) 
to a negative frequency-squared (regime II), and then back to a positive frequency-squared 
(regime III). At the interface of these regions, the frequency vanishes (w(t) = 0), and the 
condition |w/ci; 2 | <C 1 is maximally violated. In practice, this condition is violated for some 
time interval around the points where the frequency-squared changes sign. 

Assuming that the frequency is slowly varying, we begin by writing the lowest-order 
WKB approximation to the solution of Eq. 5.11 in each region [88] 


x\{t) 

Xk(t) 

x^(t) 
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V 2u} k{t) 
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\/2 uik(t) 
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A) 
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^ J uj k {t')dt'^j , 

(5.13) 
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(5.14) 
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P 

exp | 

^ J Uj k (t')dt'^j , 

(5.15) 

y/ 2uj k(t) 


where D 2 (f) = — w 2 (t), and Wfc(ti) = bJkft-i) = 0 and the integrals in the expression for 
xl/ 1 (t) are performed for the region where u 2 {t) > 0. We are ultimately interested in the 
two coefficients a and /3 which describe the amplitude of the mode after its brief growth 
due to the tachyonic instability. In regime I we match the solution to the asymptotic past 
(Bunch-Davies vacuum), giving us a® = 1 and /3q = 0. After matching at the two interfaces 
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between regimes I and II and regimes II and III (by Taylor-expanding u 2 (t) near the turning 
points and using the asymptotic form of the Airy functions) the coefficients are 

a = e Xk , and j3 = — ie Xk e~ 2ldk , (5.16) 

where 

ft 2 ft 1 

Afc = / Qk(t')dt', and 0^= ujk{t')dt'. (5-17) 

J t\ Jto 

Before proceeding, we note that the WKB approximation is accurate only in the broad reso¬ 
nance regime, g> 1. This translates to the condition that 

— %»1. (5.18) 

am j 

For typical values of the coupling and the inflaton amplitude, this condition restricts the 
validity of the WKB approximation to k > m. As shown in Fig. 5, the maximum amplifica¬ 
tion occurs well within the region of validity of our WKB calculation, especially for higher 
couplings. 

We study each polarization individually, starting with the mode A + . In the absence of 
back reaction, the axion velocity, <t>, takes its maximum positive value at time (t* — to) ~ 
4.5 m' 1 after inflation ends, 5 which sets the interval during which A + is tachyonic for the 
first time. The evolution of the inflaton field c j>[t) is close to sinusoidal at this stage, meaning 
that we can use a modified static-universe approach. However, we employ our full expanding- 
universe WKB-method as a way of testing its accuracy and providing a unified treatment of 
the two gauge-field polarizations. 

The evolution of the second polarization, A~, is more involved due to the fact that 
regime II, characterized by uj 2 < 0, starts while the Universe is still inflating and continues 
through the end of inflation into the first oscillation of the inflaton. This complication does 
not significantly alter our analysis. To proceed, we simply need to initialize the mode in 
regime I, sufficiently early during the inflationary stage, before its tachyonic transition and 
follow it, using the WKB approximation, through this transition and the end of inflation. 
The final formulas are exactly the same in this case. 

We can evaluate the validity and accuracy of the semi-analytical method described above 
by comparing the results directly with a full numerical solution of the linear equations of mo¬ 
tion. Using MATHEMATICA, we solve the homogeneous Klein-Gordon and Friedman equations 
for the evolution of the background axion and spacetime. On this background, we follow three 
gauge-field modes for At whose wavelengths are equal to the horizon, and one and two e- 
foldings smaller than the horizon at the end of inflation, i.e. k/(aH) = jl, e, e 2 } respectively 
at the end of inflation. We also track three modes for AT, one that exits the horizon one 
e-fold before the end of inflation, a mode whose wavelength is equal to the horizon, and a third 
that is one e-fold smaller than the horizon at the end of inflation, i.e. k/(aH) = {e" 1 , 1, e} 
respectively at the end of inflation. In order to facilitate the comparison, we take the mode 
amplitudes to have unit size at the start of our simulation. The WKB condition |cu/a; 2 | <C 1 
is violated around the points where u = 0 for each mode and |cU/c^ 2 1 < 0(0.1) during the 
tachyonic regime. This limits the accuracy of the approximation. The approximation could 
be improved by making use of a transformation of variables similar to [89], however, given the 

5 In this section we denote by a subscript 0 the value of a quantity at the end of the inflationary phase 
where en = —H/H 2 = 1 for the first time. 
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good agreement with the numerical results shown in Fig. 4, and the fact that we are trying 
to understand the results of full lattice simulations rather than substitute them, we do not 
attempt to refine the WKB method used here. 

In Fig. 4 we show the excellent overall agreement between the numerical results from di¬ 
rect numerical solution in MATHEMATICA and those obtained from the WKB approximation. 
As expected the approximation diverges at the points where u 2 (t) = 0, and closely follows 
the curves as one moves away from these points. Furthermore, we can see that for both A + 
and A~ the accuracy of the approximations decreases with decreasing wavenumber k (leading 
to decreasing q), especially after the tachyonic transition. Although we are pushing the limits 
of the WKB approximation, the accuracy of the resulting growth rates is remarkably robust. 





Figure 4. The amplification of gauge fields during the first tachyonic instability phase for A + (left) 
and A~ (right) based on numerical simulations and the semi-analytic calculation (for a/f = 70 mh 1 ). 
On the left panel the lines and dots correspond to numerical results for wavenumbers k/{aH) = 
l(blue), k/{aH) = e (black) and k = k/{aH ) = e 2 (red). On the right panel the lines and dots 
correspond to numerical results for wavenumbers k/(aH) = e~ 1 (blue), k/(aH) = 1 (black) and 
k/{aH) = e (red). In all cases aH is evaluated at the end of inflation, and t = 0 corresponds to the 
end of inflation. 


These results demonstrate that this semi-analytical method can be used to accurately 
estimate the growth of the gauge fields during the first inflaton oscillation, if one neglects 
rescattering and back-reaction effects. In Fig. 5, we plot the growth factor X^. which shows 
by how much the amplitude of each gauge mode has grown after the first tachyonic regime. 
Note that this WKB method breaks down at small wavenumbers. 

For the mode A + we can make a comparison to the static-universe calculation by us¬ 
ing rescaled parameters. The amplitude of the axion oscillations has decreased due to the 
expansion of the Universe after the end of inflation (as shown in Fig. 3). For the tachyonic 
half-period of interest, the behavior of the axion held is well approximated by 

cf)(t) « —0.05mmpi cos (mt + A9), (5.19) 

where A 9 is a phase offset that allows the time of zero-crossing to correspond to our model. 
The wavenumbers of the modes under consideration are also redshifted. The modes used for 
our full WKB calculation were measured in terms of the scale factor at the end of inflation 
a o, while the Universe has grown by a factor of a/a® ~ 2.6 by the middle of the first tachyonic 
regime for A + . We thus rescale the wavenumber used in our static-universe calculation by 
2.6. As we can see in Fig. 5 the results are very close, giving us a simple physical way to 
understand the result of using the WKB method in an expanding universe. 
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The polarization A~ is more complicated due to the fact that it becomes tachyonic 
during the inflationary phase. From Fig. 5 note that the growth factors for a given mode are 
significantly larger, and a larger range of wavenumbers get amplified. We do not compare this 
case with a static-universe approximation, since </> is not well approximated by a sinusoidal 
function in the regime where At is tachyonic, even after inflation has ended. This is due to 
the Hubble friction term, as shown in Fig. 3. 
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Figure 5. The growth factor X j~ versus the comoving wavenumber after the first tachyonic instability 
phase for A + (left) and A~ (right) based on the semi-analytic calculation. The dotted lines correspond 
to our WKB formulas, while the continuous lines on the left plot correspond to a modified static- 
universe approximation. The different lines correspond to different couplings a/f = 35 m^ 1 (blue 
dashed), a/f = 45 77?.“/ (red dot-dashed) and a/f = 55to“j (black dotted) for both plots. 

There is one further complication we need to keep in mind when comparing the semi- 
analytic predictions of this section with the full lattice simulations presented in Section 6. 
Once the coupling to the gauge fields is turned on, these fields begin to backreact. This causes 
inflation to end at a slightly different value of f for each coupling a/f, as shown in Table 1. 
We also show the values of the field and field velocity two e-foldings before inflation ends. Note 
that these values are all close to the zero coupling case, which indicates that backreaction 
has not yet become important yet.The WKB analysis of this section was performed using 
the evolution of the axion in the limit when we neglect back reaction from the gauge fields, 
equivalently in the limit of a/f —> 0. 


m p ia/f 

(N 

1 

O 

<^e-2 /m 

'/’end 

^end/ 

0 

0.563 

-0.158 

0.201 

-0.143 

30 

0.563 

-0.158 

0.201 

-0.143 

35 

0.563 

-0.158 

0.201 

-0.143 

40 

0.564 

-0.158 

0.203 

-0.142 

45 

0.571 

-0.158 

0.217 

-0.142 

50 

0.580 

-0.158 

0.237 

-0.147 

55 

0.579 

-0.158 

0.239 

-0.132 

60 

0.572 

-0.158 

0.235 

-0.115 

65 

0.564 

-0.158 

0.231 

-0.102 


Table 1. Field conditions two e-folds before the end of inflation and at the end of inflation for 
V = rnff 2 /2 and m/m p i = 1.06 x 10 -b including the gauge-field backreaction. Note that </> en d does 
not monotonically increase with the coupling. This is likely an artifact of the approximation. 
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6 Non-linear lattice simulations 


The results of Section 5.1 imply that at the couplings of interest, the dynamics of the axion- 
gauge field system very quickly become highly non-linear and require numerical analysis. 
Numerical methods that evolve scalar fields in an expanding background have been around 
for a couple of decades. Many numerical methods have been developed [90-93], and all have 
regimes in which they are most successful. Here we use GABE [94] (as in [59]), where the 
axion and the gauge field are defined on a discrete lattice (grid) with 256 3 points. The 
software uses a second-order Runge-Kutta integration method to solve Eqs. 2.10, 2.12 and 
2.13 alongside the self consistent expansion of space-time Eq. 2.14. We work in Lorenz gauge, 
= 0, to evolve the gauge fields. In this gauge, the Gauss law constraint becomes a 
dynamical equation for Aq which we solve in parallel with Eqn. 2.13. Note that, although 
we initialize our fields using solutions of the linear equations of motion in Coulomb gauge, 
these gauge fields are equivalent to the fields in Lorenz gauge. To see this, note that if one 
begins with fields in Coulomb gauge, and performs a gauge transformation to Lorenz gauge 
^Lorenz _ ^Coulomb _the Lorenz gauge condition (d^ H^ orenz = 0) becomes (d^—didi)a = 
0, which is simply the residual gauge symmetry of the Lorenz gauge. Unless otherwise noted, 
all simulations use a box size that is L = 15 m -1 at the end of inflation and are run using 
the parallel processing standard OpenMP with 12 threads. The equations of motion for the 
gauge degrees of freedom, An. as well as the axion, cj), are integrated using a second-order 
(midpoint method) Runge-Kutta integration scheme. We normalize quantities so that wave 
numbers are expressed in terms of the size of the box at the end of inflation, L = 15 m -1 , 
and a = 1 (at the end of inflation). 

We begin our simulations two e-foldings before inflation ends, defining this point as r = 0. 
As described at the end of Section 3, we determine the value of the homogeneous field and 
its derivative by numerically evolving (using MATHEMATICA) the system of Eqs. 3.19 - 3.21 
together with the approximations of Eq. 3.23. At this point the box-size is Lq = Le -2 ~ 2m -1 
which is just larger than the Hubble Scale, iL -1 = ^/3/87 t p -1 / 2 ~ 1.2m -1 , where the final 
approximation varies slightly from coupling to coupling. We initialize the power in the A± 
modes by numerically evolving Eq. 3.8 for each physical mode, tracking it from when it was 
well inside the horizon (r —> —oo) until two e-foldings before the end of inflation. Since we 
have no analytic form for the fluctuations of cj) at this time, we initialize it in the Bunch-Davies 
vacuum, {^k) = 1 /\/2a7, where it is an excellent approximation due to the fact that almost 
all of our modes are sub-horizon. Using this procedure, the modes that leave the horizon 
during the final two e-foldings of inflation are done so self-consistently and the spectrum of 
perturbations for 0 is consistent with our equations of motion. Fig. 6 shows the spectra at 
the beginning of the simulation and the end of inflation showing the amplification of large 
wavelength modes during the final stages of inflation. 

After we set the initial spectrum of A ^ in momentum space we project these onto the 
gauge fields 

A k = e+(k )A+ + e_(k)AjT, (6.1) 

and (inverse) Fourier transform them into configuration space using a set of projection oper¬ 
ators, efj, that satisfy the relations 

k • e± = 0, k x e± = (6.2) 

These relations set only the spatial components of the gauge field, A(x, r = 0), on the initial 
surface. Since we are numerically tracking the values of the full four-potential, we must 
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Figure 6. Upper panels: The initial power spectrum of the fields, A + (red, solid) A- (blue, dashed) 
and 4> (black, dotted) at a/f = 45 m^ 1 . The left panel shows the spectra at the beginning of the 
simulation which occurs at two-efoldings before the end of inflation, a = e~ 2 . The right panel shows 
the spectra of the same three fields at the end of inflation. The lower panels show the same spectra 
at a coupling of a/f = 65 771“^, the highest coupling that we simulate. Note the long-wavelength 
fluctuations of the axion generated during the end of inflation by for the largest couplings. 


check to make sure the Lorenz gauge condition, = 0, is obeyed in configuration space 

as required by our equations of motion. The definition of the polarizations, Eq. 6.2, requires 
that A® = 0 (A 0 = constant) on the initial slice. Therefore any choice of A± (with the choice 
h4o = 0) obeys the gauge condition. We are, of course, neglecting any effect that the initial 
conditions of A^ have on </> on the initial surface. However, we begin our simulations during 
inflation where all modes of interest are sub-horizon. After the fields are initialized, there 
is an amplification of modes of <ft that reaches equilibrium well before the first zero-crossing 
of the held. To make sure that we don’t depart from the gauge-condition surface, obeying 
Lorenz gauge, we track the size of 


G(r) 


_ d 0 A 0 — V • A _ 

^/(d 0 A 0 y + (9iAr)2 + (5 2 A 2 )2 + (5 3 A 3 ) 2 


(6.3) 


as in [59]. Fig. 7 shows a plot of G(r) averaged over the box and the RMS averaged value, 
yf (G 2 ), for a simulation with a/f = 45 m^ 1 . We compare this at two resolutions, N = 128 
and N = 256. The reader should understand this as a proxy for staying on the gauge surface— 
when the simulation diverges from satisfying the gauge condition, this measure quickly grows 
large. At early times, since we do not ‘cutoff’ the initial spectra, the power in the very highest 
frequency modes contribute to variations in the finite-derivatives we use to calculate G(r) are 
artificially important, and contribute to the size of the RMS value, since the held values are so 


19 















small. Once the modes of the gauge held are populated, at the end of inflation, this measure 
is increasingly good. 




Figure 7. The left panel shows the four terms of that contribute to calculating the gauge 

condition at an arbitrary point in the simulation box: 8qA° (black), d\A 1 (red), c^A 2 (blue), and 
c? 3 d 3 (purple). The right panel shows a plot of G(t) averaged over all points on the box, (G), (black, 
solid) and the RMS average, \J (G 2 ), (red, solid), at N = 256 compared to the same two quantities 
(same colors but dashed) for N = 128. This parameterizes how well the gauge condition is satisfied 
and is shown for a/f = 45 m~ j 

The goals of our simulations are to see if the tachyonic and/or parametric instabilities 
identified in Section 5 lead to the efficient generation of gauge modes, whether the energy 
deposited in those gauge modes is enough to preheat the Universe, and whether that final state 
has any anisotropy in the two polarization states of the gauge held, as a naive interpretation 
of Fig. 5 would suggest. To parameterize the success of the first and second of these questions, 
we calculate the total energy in the gauge held (as defined in Eq. 3.20) 

Pem = \{ e2 + B 2 ) , (6-4) 

although we express this as a fraction of the total energy density of the Universe, Pem/ ptot, 
for various couplings. 

As the simulations progress, we see that energy is, generically speaking, quickly and 
efficiently transferred into the gauge helds. Fig. 8 shows the evolution of the ratio of the 
energy in the gauge helds to the total energy as a function of time through the simulation for 
different values of the couplings. In all cases where a/f > the Universe completely 

preheats and almost all of the energy of the simulation is transferred into the gauge helds. 
In most cases, this happens during the hrst few oscillations of the axion held, justifying the 
need for full non-linear simulations. Only in the case of a marginal coupling a/f = 45 m^ 1 
does the axion take almost two full oscillations before the gauge held dominates the energy 
density. Below, we discuss a few regimes identihed by these simulations. 

6.1 Early tachyonic resonance 

The most significant difference between this model and previous studies of gauge helds during 
reheating following inflation [59] is the prediction that one of the polarizations of the gauge 
held is tachyonic each time the homogeneous mode crosses zero. This should persist as long 
as the inhaton is (dominantly) coherent and we can treat it as a strictly time-dependent 
quantity in Eq. 5.3. We are using a box whose physical size increases over the course of the 
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Figure 8. The total energy fraction in the gauge field pEM/ptotai as a function of time for a variety 
of couplings. We probe couplings between a/f = 35 m^ 1 to a/f = 65 to^ 1 in increments of a/f = 
JjnTpi 1 - The left panel shows these couplings as a function of time for each run. The couplings go 
from largest on the top (red) to the lowest on the bottom (blue) along the rainbow spectrum. The 
two lowest couplings, a/f = 35and 40 m~/ do not completely preheat. The right panel shows 
the maximum ratio of /?EM/ptotai for each value of the coupling, a/f. The initial energy of the gauge 
fields is red-shifted away during the last two e-foldings of inflation since we are unable to introduce 
shorter-wavelength modes during the simulation. 


simulation. The longest wavelength we are able to probe at the beginning of the simulation 
corresponds to a minimum wavenumber, fc m ; n = 27r/L ~ 0.4 mT 1 (although this does grow 
during the simulation as the scale factor increases). 

During the first oscillation, we expect low frequency modes of the A~ mode to be excited. 
This should be extremely efficient up to a maximum wavenumber, (k/a) ~ m(j>o(a/ /). In 
each simulation we have a slightly different value of the inffaton field at the end of inflation, 
<f> o, depending on the value of the coupling a/f - although they vary only by only about 
15%, as shown in Table 1. Fig. 9 is a comparison of the strength of this tachyonic instability 
including early effects of backreaction and rescattering (which seems to be largely missing 
from analytic studies of preheating in these models) showing the effect of this first tachyonic 
regime on the power in the gauge field. 


cn_ 

I 




Figure 9. We compare the power spectrum of the A~ mode for simulations in which the coupling is 
a/f = 35 771^ (red, solid) a/f = 45 m~/ (blue, dashed), and a/f = 55171/,/ (black, dotted). The left 
panel shows the power spectrum of the A~~ at the time when the axion crosses zero for the first time 
for these three couplings. The right panel shows the ratio of A~ evaluated at the first zero crossing to 
the initial spectrum, capturing the enhancement of A~ from the beginning of the simulation to this 
time. 
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For moderate couplings, a/f > 50 m^ 1 we see that this first wave of tachyonic instability 
is extremely efficient and an 0(1) fraction of the total energy density of the Universe is 
already deposited into the gauge fields. This is a highly-tachyonic regime. By the time 
that the homogeneous mode of the axion has completed its second zero crossing, it is already 
incoherent and the majority of the energy of the axion is in higher-momentum modes. Fig. 10 
shows the final spectra of the gauge field for one of these cases. Note that the spectra do 
not exhibit any discernible signs of band structure. This is expected from the analysis of 
Section 5.1. During the first tachyonic regime all modes within the tachyonic window get 
amplified in the first oscillation. The band structures form due to resonances after multiple 
oscillations. However, in these cases, multiple coherent oscillations of the axion background 
are prevented due to the strong gauge-field back reaction. This back reaction is so strong 
that at the largest couplings probed here the homogeneous axion condensate does not cross 
zero. We call this period the early tachyonic regime; the period of strong tachyonic growth 
that occurs during the first oscillation of the axion. 

Naively, one might have predicted that these larger couplings would have lead to a 
preheated state that is highly polarized. The lack of oscillations of the axion effectively 
prevents the tachyonic regime of the A + mode from developing. On the other hand, the 
A~ mode is strongly produced as the axion condensate relaxes to zero. Further, once the 
axion condensate has become unimportant, any energy transfer essentially ceases. While this 
is certainly the case for the process of interest, during preheating, the excited A~ modes 
strongly re-scatter off (j) and source A + modes. This process is very efficient on sub-horizon 
scales, and the resulting spectra show very little difference in power between each polarization 
on these scales. This effect is demonstrated in Fig. 10, where we plot the final spectra of the 
gauge modes, A±, and the axion, (j), in the full model. We also show the effect of artificially 
switching-off this rescattering by eliminating the terms in the equation of motion for the axion 
that couple it to the gauge field, the term proportional to a/f in Eq. 2.10. This prevents non- 
linearities from developing in the axion and the subsequent rescattering of the gauge modes. 
Note that the final state in these cases is strongly polarized, which demonstrates the efficiency 
of re-scattering. In the case where we artificially block this back-reaction, our simulations do 
not conserve energy and the Hubble parameter rises during the tachyonic phases; this causes 
the spectra for the gauge fields to be larger than expected in the full simulations by a small 
factor. The right panel of Fig. 10 (as well as Figs. 11 and 12) should be considered illustrative. 
Comparing this to Fig. 5 we see that in the absence of back-reaction the amplitudes of two 
polarizations differ by several orders of magnitude, as expected by our semi-analytic study. 

Even in the presence of back-reaction, the spectra are still polarized in the infrared (the 
spectra here differ by a couple orders of magnitude, see left panel of Fig. 10), suggesting that 
the asymmetric long-wavelength modes might have some observational consequence. 

When the tachyonic growth factor is very large (Fig. 5), the energy transferred to the 
gauge held (A~) after the first oscillation is comparable to the initial energy stored in the 
inhaton condensate, and the Universe preheats almost instantaneously. The couplings that 
lead to instantaneous preheating can be estimated from the results of Section 5.1. The energy 
density of the background inhaton is 

p<t> = * (dtcfi) 2 + 1m 2 q i> 2 ~ 0.05 m 2 ?^! (6.5) 

where we used the value of 4> a t the end of the hrst tachyonic stage. The energy density of 
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Figure 10. The final power spectra of A + (red, solid), A (blue, dashed) and </> (black, dotted) for a 
simulation in which the coupling is a/f = 60771“/. Even though the tachyonic process is only efficient 
for one of the polarizations, the A ~, of the gauge field, there is significant back-reaction onto the 
modes of (j> which source modes of A + . The left panel shows the final spectrum for a simulation of our 
model. The right panel shows the same simulation, however, with the interaction term eliminated from 
the equations of motion for the axion field. This comparison shows that a lot of power is transferred 
between polarizations, mediated by the axion. Both panels are evaluated at a = 8. 


the gauge field can be calculated from Eq. 6.4 to be 


Pem = ~j d 3 k(\d T A \ 2 + k 2 \A | 2 ) 
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where we dropped the time-derivative term, since the mode functions do not grow between 
tachyonic regions. Further we only considered the A~~ mode, since during the first tachyonic 
region the A + mode is exponentially smaller, hence it can be safely ignored in the calculation 
of the energy density. The last approximate equality is based on the fact that the function 
e x k i s sharply peaked. Using the values given in Fig. 5, we get Pem/ P<f> ~ 10“ 4 ,10“ 2 ,10 2 
for mpia// = 35,45,55 respectively after the first tachyonic regime at a(t ) ~ 2. The first 
two values agree with Fig. 8. The analytical result for mp\ot/f = 55 signifies that we have 
entered the region of strong energy transfer and strong back-reaction, hence we cannot treat 
the axion field as decoupled. These estimations immediately show the region of couplings, for 
which instantaneous preheating can occur, where full non-linear simulations are unavoidable. 

There are additionally cases, however, in which the tachyonic instabilities are present, 
but not sufficiently efficient to deposit an 0( 1) fraction of the total energy into the gauge 
field during the first oscillation. In these cases, it can take up to ten oscillations of the 
homogeneous mode of 4> before the Universe is preheated. It is still possible for the Universe 
to totally preheat in these cases, it just takes longer. This regime also results in an unpolarized 
final state, as seen in Fig. 11, where the axion, again, plays a significant role in balancing the 
two helicity states. 


6.2 Parametric resonance 

For lower values of the coupling, a/f, the efficiency of the early tachyonic regime is not high 
enough to completely preheat the Universe. At the end of inflation, in the homogeneous field 
limit, the potential and kinetic energy of the inflaton are approximately equal and so 
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Figure 11. The final power spectra of A + (red, solid), A~ (blue, dashed) and f (black, dotted) for 
a simulation in which the coupling is a/f = 45 m ~/. In this case, the final spectrum is not polarized; 
the power in each of the two modes is almost identical. The left panel shows the final spectrum for a 
simulation of our model. The right panel shows the same simulation, however, with the interaction 
term eliminated from the equations of motion for the axion field. This comparison shows that the 
transfer of power between modes is still important in this regime. Both panels are evaluated at a = 8. 


which defines the smallest wavelength that is dynamical. Since the maximum wavenumber 
that can be amplified is 

*77? = (6-8) 

a(t) f 

we see that the band that gets amplified shrinks as a/f gets smaller. Examining Fig. 5, we 
see that both the growth factor, as well as the regime of amplified wavenumbers, shrink as 
the coupling a/f gets smaller and the early tachyonic regime is not sufficient to transfer most 
of the inflaton energy into the gauge fields. 

While at low couplings, gauge-held production during the early tachyonic-regime is not 
strong enouigh to reheat the Universe, the persistence of coherent oscillations of the axion 
condensate allows for parametric resonance. Parametric resonance continues to deposit power 
into the gauge held - independently of polarization- for many oscillations. In Fig. 12 we show 
the power spectra of the gauge held polarizations when a/f = 30 m^ 1 . Since many different 
modes have passed through different parametric-resonance bands over many oscillations, a 
complex spectral structure forms. However, this spectral structure is shared between the 
two polarizations equally, because in the limit of many oscillations the two polarizations 
obey identical equations. While the final state is polarized, this is a residual effect of the 
initial conditions. The formation of these spectral bands can be explained in the WKB 
approximation (for g > 1) by a process analogous to multiple scatterings from a periodic 
potential leading to constructive and destructive interference (as discussed in [76]), as well as 
through Floquet theory (for q < 1). 

However, the question that arises by inspecting Fig. 8 is the inability of lower couplings 
( a/f < dOmpj 1 ) to fully reheat the Universe. The parameters of the Mathieu equation of Eq. 
5.6 in an expanding matter-dominated universe become 


( k V I 2 _ fc q (fo 

\ma(t)) 91 2 ’ ^ ma(t) f t 


(6.9) 


where we took the scale factor to evolve as a(t) ~ t 2 / 3 and the envelope of the inflaton 
oscillations to decay as (f>o(t) = </>oa~ 3 / 2 = (f>o i _1 - We thus see that as time progresses 
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the product of the coupling with the axion amplitude becomes smaller. Since this product 
controls the strength of the parametric resonance, we can call it the “effective coupling”. For 
a/f = dOmpj 1 , the effective coupling becomes unity at t ~ 9 m or equivalently a(t) ~ 4. The 
analysis of [38] shows that a pseudo-scalar inflaton derivatively coupled to U( 1) gauge fields 
with 0(1) coupling does not fully preheat, since the growth rate is similar in strength to the 
red-shifting of the gauge field amplitude due to the expansion of the Universe. Transferring 
this known result to our case simply means that, in models where most of the inflaton energy 
has not been transferred to the gauge fields by the time the effective coupling equals unity, 
no significant preheating occurs. This both explains the behavior seen in Fig. 8 as well as 
further demonstrates the importance of the first few axion oscillations for preheating into 
gauge fields. In any case, perturbative decay of the axion to gauge field is still present and 
can eventually reheat the Universe, as discussed in Section 4. 



Figure 12. The final power spectra of A + (red, solid), A~ (blue, dashed) and 4> (black, dotted) for 
a simulation in which the coupling is a/f = 35 . In this case, we see a more complex spectrum 

whose modes have been amplified by a number of (varying) resonance bands. The left panel shows the 
final spectrum for a simulation of our model. The right panel shows the same simulation, however, 
with the interaction term eliminated from the equations of motion for the axion field. In both panels 
the initial asymmetry between the power in the two helicity modes persists until the end of the 
simulation. Both panels are evaluated at a = 8. 


6.3 Monodromy potential 

It is worth verifying that the results presented here are (at least somewhat) insensitive to the 
specific form of our axion potential. To test this, we study reheating following monodromy 
inflation in the potential of Eq. 2.3. This potential introduces another scale to the problem 
and causes the oscillations to be increasingly anharmonic when the field probes the region 
4> > (f> c . As noted above, the fluctuations observed in the CMB are insensitive to the value of 
(f> c (as in Section 2) and thus we can treat it as a free parameter. For this section we choose 
4> c = 0.02 m p p We then chose a set of couplings, a/ f, and checked to see how much energy is 
transferred into the coupled gauge field. Fig. 13 shows the effect ou the reheating efficiency 
for this potential. 

As can be seen from Fig. 13, monodromy inflation is less efficient at reheating into 
gauge fields than quadratic, ?n 2 0 2 , inflation, typically requiring larger couplings to achieve 
reheating within the first axion oscillation. This is not particularly surprising, as the axion 
typically rolls at a slower rate relative to the Hubble rate on this potential, and consequently 
gauge-field production is lower. This is evidenced by the larger allowed values of the coupling 
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a// by the black-hole abundance bounds, which are allowed to be around approximately 50% 
larger than the m?(j) 2 couplings. We find similar behavior in the reheating epoch, the lowest 
coupling that completely reheats in this case is a/f ~ 70 m~^ which should be compared with 
a/7 45m.pi 1 f° r 777,2 0 2 - Tlhs result is due to the fact that the efficiency of the tachyonic 

regime is sensitive to the axion velocity. 



Figure 13. The left panel shows the total energy fraction in the gauge field Pem/ptotai as a function 
of time for a variety of couplings. These range from a/f = 120 m^ 1 to a/f = 30 to^ 1 in increments 
of a/f = lOjn.pj 1 . The couplings go from largest on the top (red) to the lowest on the bottom (blue) 
along the rainbow spectrum. The right panel shows the fraction of the axion energy contained in 
regions where p$ > 4 (p^,) (defined by Eq. 6.10) using the same color scheme. This is a proxy to test 
for the relative abundance of oscillons in these simulations. 


It is known that under certain conditions pseudo-stable inflaton lumps, or oscillons 
[81, 95-97], can form from the post-inflationary detritus. The interplay between scalar and 
gauge fields in oscillons has only been studied for models where the scalar held is charged 
and the gauge held is non-Abelian [98, 99]. The oscillons that form from the fragmentation 
of the inhaton following inflation can be long lived and have drastic consequences for theories 
of the early Universe. This is because they behave like pressure-less dust and can change 
the early expansion history. Our methodology is sufficiently robust to allow us to study the 
creation and decay of these structures, as well as to probe the effect of these structures on our 
reheating history. We test for the formation of oscillons in our simulations, by computing the 
fraction of energy contained in regions of high density compared to the average background 
energy density. Specifically, we compute the total energy in regions where the local density 
is greater than four times the average density, and then compare this to the total energy in 
the box [81], 


fp*> 4<P 0 > P* dV 

f P<p dV 


( 6 . 10 ) 


We note that our definition differs from that in [81], since we only calculate the energy in 
the axion held, when we compute the numerator and denominator and we restrict the 
integral in the numerator to regions of very high energy density. We also only integrate the 
numerator for regions where the axion energy density is four times the average axion energy 
density (instead of twice the axion energy density). The right hand panel of Fig. 13 shows 
this statistic and shows that the time at which oscillons form is highly dependent on the 
coupling, suggesting that the gauge fields play a role in the creation (and possible decay) of 
these structures. We delay a full treatment of oscillons for a future publication. 
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7 Conclusions 


In this work we have studied preheating into 17(1) gauge fields following axion inflation. 
The shift-symmetry of the axion means that the form of the couplings to matter fields is 
highly restricted. A class of allowed couplings is via derivative interactions, and consequently 
dimension-five couplings to gauge fields via Chern-Simons terms are expected in the theory 
from an effective held theory point of view. The size of such a coupling is unknown, however, 
the effect of rescattering of helically polarized gauge bosons off the axion condensate and the 
subsequent generation of curvature fluctuations during the inhationary epoch places an upper 
bound. In this work we have explored the phenomenological consequences of such a coupling 
on the post-inhationary evolution of axion-driven inhationary models. 

Using lattice simulations and semi-analytic methods we have shown that preheating 
into Abelian gauge helds via these Chern-Simons interactions can be extremely efficient. 
In particular, at the middle to upper range of the couplings allowed by constraints due to 
over-production of primordial black holes, we find that reheating is essentially instantaneous, 
proceeding via a phase of early tachyonic resonance and completing within a single oscillation 
of the axion. The resulting Universe ends up in an un-polarized state due to strong rescattering 
effects on scales that are sub-horizon during reheating. Scattering of amplified gauge held 
modes into axion fluctuations generates the second polarization extremely efficiently. On 
super-horizon scales, the asymmetry that develops due to the tachyonic instability of one of the 
gauge modes during inflation, remains. The Universe in these high coupling cases is radiation 
dominated and is characterized by very efficient preheating where non-linear dynamics and 
back reaction become important almost immediately. As the coupling is decreased, this 
phase of early tachyonic resonance is weakened, and the axion oscillates multiple times before 
preheating completes. During these multiple oscillations equal levels of both polarizations 
of the gauge held are excited due to tachyonic and parametric resonance. Decreasing the 
coupling further yields a brief window where parametric-resonance effects become important, 
before preheating abruptly shuts off and non-linear effects cease to be important. At these 
lower couplings, non-linear effects are never important and the Universe reheats perturbatively 
due to the decay of the axion into gauge bosons with a reheating temperature near 10 9 GeV. 
Note that all of the couplings we have considered in this work are far below those values that 
give any observable effects in the CMB. At the couplings that saturate the bounds considered 
in the recent Planck paper [100], reheating would be essentially instantaneous and, perhaps, 
accompanied by an overproduction of primordial black holes. 

We studied reheating in two different axion potentials: the simplest model of chaotic 
inflation with a monomial potential, m 2 ^ 2 , and the simplest model of axion monodromy infla¬ 
tion. Phenomenologically, these two potentials have similar (p)reheating behavior. However, 
generically for a given coupling, the monodromy potential has a lower reheating efficiency 
compared to the chaotic inflation case due to the slower initial axion velocity. The anhar- 
monic nature of the axion increases the features of preheating in this model and warrants 
further study. 

For the case of monodromy inflation, it is well known that at the end of inflation pseudo¬ 
stable classical lumps of the axion held - oscillons - can form and lead to a period of matter 
domination before the onset of reheating. Our numerical investigations suggest that these 
couplings to gauge helds only strengthen the formation of these oscillons and leads to their 
formation at an earlier epoch compared to the uncoupled case. However, the Universe remains 
radiation dominated due to the bath of gauge bosons produced. We leave a full study of 
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oscillons in these models to future work. 


Acknowledgments 

We thank Mustafa Amin for useful discussions. PA gratefully acknowledges support from 
a Starting Grant of the European Research Council (ERC STG grant 279617) and the hos¬ 
pitality of the Department of Applied Mathematics and Theoretical Physics at the Univer¬ 
sity of Cambridge where some of this work was completed. JTG and TRS are supported 
by the National Science Foundation, PHY-1414479. We acknowledge the National Science 
Foundation, the Research Corporation for Science Advancement and the Kenyon College De¬ 
partment of Physics for providing the hardware used to carry out these simulations. EIS 
gratefully acknowledges support from a Fortner Fellowship at the University of Illinois at 
Urbana-Champaign. 

References 

[1] BICEP2 Collaboration Collaboration, P. Ade et al., Detection of B-Mode Polarization at 
Degree Angular Scales by BICEP2, Phys.Rev.Lett. 112 (2014) 241101, [arXiv: 1403.3985]. 

[2] Planck Collaboration Collaboration, R. Adam et al., Planck intermediate results. XXX. 

The angular power spectrum of polarized dust emission at intermediate and high Galactic 
latitudes, arXiv: 1409.5738. 

[3] BICEP2 Collaboration, Planck Collaboration Collaboration, P. Ade et al., A Joint 
Analysis of BICEP2/Keck Array and Planck Data, Phys.Rev.Lett. (2015) 

[arXiv:1502.00612], 

[4] V. F. Mukhanov and G. V. Chibisov, Quantum Fluctuation and Nonsingular Universe. (In 
Russian), JETP Lett. 33 (1981) 532-535. 

[5] D. H. Lyth, What would we learn by detecting a gravitational wave signal in the cosmic 
microwave background anisotropy?, Phys.Rev.Lett. 78 (1997) 1861-1863, [hep-ph/9606387]. 

[6] K. Freese, J. A. Frieman, and A. V. Olinto, Natural inflation with pseudo - Nambu-Goldstone 
bosons, Phys.Rev.Lett. 65 (1990) 3233-3236. 

[7] F. C. Adams, J. R. Bond, K. Freese, J. A. Frieman, and A. V. Olinto, Natural inflation: 
Particle physics models, power law spectra for large scale structure, and constraints from 
COBE, Phys. Rev. D47 (1993) 426-455, [hep-ph/9207245]. 

[8] T. Banks, M. Dine, P. J. Fox, and E. Gorbatov, On the possibility of large axion decay 
constants, JCAP 0306 (2003) 001, [hep-th/0303252]. 

[9] S. Dimopoulos, S. Kachru, J. McGreevy, and J. G. Wacker, N-flation, JCAP 0808 (2008) 003, 
[hep-th/0507205], 

[10] R. Easther and L. McAllister, Random matrices and the spectrum of N-flation, JCAP 0605 
(2006) 018, [hep-th/0512102] . 

[11] T. C. Bachlechner, M. Dias, J. Frazer, and L. McAllister, A New Angle on Chaotic Inflation, 
arXiv:1404.7496. 

[12] A. R. Liddle, A. Mazumdar, and F. E. Schunck, Assisted inflation, Phys.Rev. D58 (1998) 
061301, [astro-ph/9804177]. 

[13] J. E. Kim, H. P. Nilles, and M. Peloso, Completing natural inflation, JCAP 0501 (2005) 005, 
[hep-ph/0409138]. 


- 28 - 


[14] C. Long, L. McAllister, and P. McGuirk, Aligned Natural Inflation in String Theory, 

Phys.Rev. D90 (2014) 023501, [arXiv: 1404.7852], 

[15] C. Burgess and D. Roest, Inflation by Alignment, arXiv: 1412.1614. 

[16] L. McAllister, E. Silverstein, and A. Westphal, Gravity Waves and Linear Inflation from 
Axion Monodromy, Phys.Rev. D82 (2010) 046003, [arXiv: 0808. 0706]. 

[17] E. Silverstein and A. Westphal, Monodromy in the CMB: Gravity Waves and String Inflation, 
Phys.Rev. D78 (2008) 106003, [arXiv: 0803. 3085]. 

[18] L. McAllister, E. Silverstein, A. Westphal, and T. Wrase, The Powers of Monodromy, JHEP 
1409 (2014) 123, [arXiv: 1405.3652], 

[19] N. Kaloper, A. Lawrence, and L. Sorbo, An Ignoble Approach to Large Field Inflation, JCAP 
1103 (2011) 023, [arXiv: 1101. 0026], 

[20] F. Marchesano, G. Shiu, and A. M. Uranga, F-term Axion Monodromy Inflation, JHEP 1409 
(2014) 184, [arXiv: 1404.3040], 

[21] R. Blumenhagen and E. Plauschinn, Towards Universal Axion Inflation and Reheating in 
String Theory, Phys.Lett. B736 (2014) 482-487, [arXiv: 1404.3542]. 

[22] A. Hebecker, S. C. Kraus, and L. T. Witkowski, D7-Brane Chaotic Inflation, Phys.Lett. B737 
(2014) 16-22, [arXiv: 1404.3711], 

[23] Y.-F. Cai, F. Chen, E. G. M. Ferreira, and J. Quintin, A new model of axion monodromy 
inflation and its cosmological implications, arXiv: 1412.4298. 

[24] E. Pajer and M. Peloso, A review of Axion Inflation in the era of Planck, Class. Quant. Grav. 
30 (2013) 214002, [arXiv:1305.3557], 

[25] D. Baumann and L. McAllister, Inflation and String Theory, arXiv: 1404 . 2601. 

[26] M. A. Amin, M. P. Hertzberg, D. I. Kaiser, and J. Karouby, Nonperturbative Dynamics Of 
Reheating After Inflation: A Review, arXiv: 1410.3808. 

[27] S. M. Carroll and G. B. Field, The Einstein equivalence principle and the polarization of radio 
galaxies, Phys.Rev. D43 (1991) 3789. 

[28] W. D. Garretson, G. B. Field, and S. M. Carroll, Primordial magnetic fields from 
pseudoGoldstone bosons, Phys.Rev. D46 (1992) 5346-5351, [hep-ph/9209238]. 

[29] T. Prokopec, Cosmological magnetic fields from photon coupling to fermions and bosons in 
inflation, astro-ph/0106247. 

[30] M. M. Anber and L. Sorbo, Naturally inflating on steep potentials through electromagnetic 
dissipation, Phys.Rev. D81 (2010) 043534, [arXiv: 0908. 4089]. 

[31] N. Barnaby, R. Namba, and M. Peloso, Phenomenology of a Pseudo-Scalar Inflaton: Naturally 
Large Nongaussianity, JCAP 1104 (2011) 009, [arXiv: 1102.4333]. 

[32] N. Barnaby, E. Pajer, and M. Peloso, Gauge Field Production in Axion Inflation: 
Consequences for Monodromy, non-Gaussianity in the CMB, and Gravitational Waves at 
Interferometers, Phys.Rev. D85 (2012) 023525, [arXiv: 1110.3327]. 

[33] R. Z. Ferreira and M. S. Sloth, Universal Constraints on Axions from Inflation, JHEP 1412 
(2014) 139, [arXiv: 1409.5799], 

[34] P. Adshead, E. Martinec, and M. Wyman, Gauge fields and inflation: Chiral gravitational 
waves, fluctuations, and the Lyth bound, Phys. Rev. D88 (2013), no. 2 021302, 

[arXiv: 1301 .2598], 

[35] M. Shiraishi, A. Ricciardone, and S. Saga, Parity violation in the CMB bispectrum by a 
rolling pseudoscalar, JCAP 1311 (2013) 051, [arXiv: 1308.6769]. 


- 29 


[36] J. L. Cook and L. Sorbo, An inflationary model with small scalar and large tensor 
nongaussianities, JCAP 1311 (2013) 047, [arXiv: 1307.7077]. 

[37] R. H. Brandenberger, A. Knauf, and L. C. Lorenz, Reheating in a Brane Monodromy Inflation 
Model, JHEP 0810 (2008) 110, [arXiv: 0808. 3936], 

[38] C. Armendariz-Picon, M. Trodden, and E. J. West, Preheating in derivatively-coupled 
inflation models, JCAP 0804 (2008) 036, [arXiv:0707.2177]. 

[39] J. Braden, L. Kofman, and N. Barnaby, Reheating the Universe After Multi-Field Inflation, 
JCAP 1007 (2010) 016, [arXiv: 1005. 2196], 

[40] N. Barnaby and M. Peloso, Large Nongaussianity in Axion Inflation, Phys.Rev.Lett. 106 
(2011) 181301, [arXiv: 101 1.1500], 

[41] A. Linde, S. Mooij, and E. Pajer, Gauge field production in supergravity inflation: Local 
non-Gaussianity and primordial black holes, Phys.Rev. D87 (2013), no. 10 103506, 

[arXiv: 1212. 1693], 

[42] E. Bugaev and P. Klimai, Axion inflation with gauge field production and primordial black 
holes, Phys.Rev. D90 (2014), no. 10 103501, [arXiv: 1312.7435]. 

[43] J. H. Traschen and R. H. Brandenberger, Particle Production During Out-of-equilibrium 
Phase Transitions, Phys.Rev. D42 (1990) 2491-2504. 

[44] L. Kofman, A. D. Linde, and A. A. Starobinsky, Reheating after inflation, Phys.Rev.Lett. 73 
(1994) 3195-3198, [hep-th/9405187], 

[45] J. Garcia-Bellido and A. D. Linde, Preheating in hybrid inflation, Phys.Rev. D57 (1998) 
6075-6088, [hep-ph/9711360], 

[46] S. Khlebnikov and I. Tkachev, Relic gravitational waves produced after preheating , Phys.Rev. 
D56 (1997) 653-660, [hep-ph/9701423], 

[47] B. R. Greene, T. Prokopec, and T. G. Roos, Inflaton decay and heavy particle production with 
negative coupling, Phys.Rev. D56 (1997) 6484-6507, [hep-ph/9705357]. 

[48] M. Parry and R. Easther, Preheating and the Einstein field equations, Phys.Rev. D59 (1999) 
061301, [hep-ph/9809574], 

[49] B. A. Bassett, D. I. Kaiser, and R. Maartens, General relativistic preheating after inflation, 
Phys.Lett. B455 (1999) 84-89, [hep-ph/9808404]. 

[50] J. Garcia-Bellido, Preheating the universe in hybrid inflation, hep-ph/9804205. 

[51] R. Easther and M. Parry, Gravity, parametric resonance and chaotic inflation, Phys.Rev. D62 
(2000) 103503, [hep-ph/9910441]. 

[52] A. R. Liddle, D. H. Lyth, K. A. Malik, and D. Wands, Superhorizon perturbations and 
preheating, Phys.Rev. D61 (2000) 103509, [hep-ph/9912473]. 

[53] F. Finelli and S. Khlebnikov, Metric perturbations at reheating: The Use of spherical 
symmetry, Phys.Rev. D65 (2002) 043505, [hep-ph/0107143]. 

[54] B. A. Bassett, S. Tsujikawa, and D. Wands, Inflation dynamics and reheating, Rev.Mod.Phys. 
78 (2006) 537-589, [astro-ph/0507632], 

[55] D. I. Podolsky, G. N. Felder, L. Kofman, and M. Peloso, Equation of state and beginning of 
thermalization after preheating, Phys.Rev. D73 (2006) 023501, [hep-ph/0507096]. 

[56] H. L. Child, J. Giblin, John T., R. H. Ribeiro, and D. Seery, Preheating with Non-Minimal 
Kinetic Terms, Phys.Rev.Lett. Ill (2013) 051301, [arXiv: 1305. 0561]. 

[57] A. Rajantie and E. J. Copeland, Phase transitions from preheating in gauge theories, 

Phys.Rev.Lett. 85 (2000) 916, [hep-ph/0003025]. 


30 - 


[58] E. J. Copeland, S. Pascoli, and A. Rajantie, Dynamics of tachyonic preheating after hybrid 
inflation, Phys.Rev. D65 (2002) 103517, [hep-ph/020203l]. 

[59] J. T. Deskins, J. T. Giblin, and R. R. Caldwell, Gauge Field Preheating at the End of 
Inflation, Phys.Rev. D88 (2013), no. 6 063530, [arXiv: 1305. 7226]. 

[60] J. Garcia-Bellido, D. Y. Grigoriev, A. Kusenko, and M. E. Shaposhnikov, Nonequilibrium 
electroweak baryogenesis from preheating after inflation, Phys.Rev. D60 (1999) 123504, 

[hep -ph/9902449]. 

[61] J. Smit and A. Tranberg, Chern-Simons number asymmetry from CP violation at electroweak 
tachyonic preheating, JHEP 0212 (2002) 020, [hep-ph/0211243]. 

[62] A. Tranberg and J. Smit, Baryon asymmetry from electroweak tachyonic preheating, JHEP 
0311 (2003) 016, [hep-ph/0310342]. 

[63] J.-I. Skullerud, J. Smit, and A. Tranberg, W and Higgs particle distributions during 
electroweak tachyonic preheating, JHEP 0308 (2003) 045, [hep-ph/0307094]. 

[64] J. Garcia-Bellido, M. Garcia-Perez, and A. Gonzalez-Arroyo, Chern-Simons production during 
preheating in hybrid inflation models, Phys.Rev. D69 (2004) 023504, [hep-ph/0304285]. 

[65] B. van Tent, J. Smit, and A. Tranberg, Electroweak scale inflation, inflaton Higgs mixing and 
the scalar spectral index, JCAP 0407 (2004) 003, [hep-ph/0404128]. 

[66] M. van der Meulen, D. Sexty, J. Smit, and A. Tranberg, Chern-Simons and winding number 
in a tachyonic electroweak transition, JHEP 0602 (2006) 029, [hep-ph/0511080]. 

[67] A. Tranberg, J. Smit, and M. Hindmarsh, Simulations of cold electroweak baryogenesis: Finite 
time quenches, JHEP 0701 (2007) 034, [hep-ph/0610096]. 

[68] J.-F. Dufaux, D. G. Figueroa, and J. Garcia-Bellido, Gravitational Waves from Abelian Gauge 
Fields and Cosmic Strings at Preheating, Phys.Rev. D82 (2010) 083518, [arXiv: 1006. 0217]. 

[69] J. Garcia-Bellido, D. G. Figueroa, and J. Rubio, Preheating in the Standard Model with the 
Higgs-Inflaton coupled to gravity, Phys.Rev. D79 (2009) 063531, [arXiv:0812.4624]. 

[70] A. Diaz-Gil, J. Garcia-Bellido, M. Garcia Perez, and A. Gonzalez-Arroyo, Magnetic field 
production during preheating at the electroweak scale, Phys.Rev.Lett. 100 (2008) 241301, 
[arXiv: 0712. 4263], 

[71] A. Diaz-Gil, J. Garcia-Bellido, M. G. Perez, and A. Gonzalez-Arroyo, Primordial magnetic 
fields from preheating at the electroweak scale, JHEP 0807 (2008) 043, [arXiv: 0805. 4159]. 

[72] A. Mazumdar and H. Stoica, Exciting gauge field and gravitons in a brane-anti-brane 
annihilation, Phys.Rev.Lett. 102 (2009) 091601, [arXiv: 0807. 2570]. 

[73] R. Allahverdi, A. Ferrantelli, J. Garcia-Bellido, and A. Mazumdar, Non-perturbative 
production of matter and rapid thermalization after MSSM inflation, Phys.Rev. D83 (2011) 
123507, [arXiv: 1103.2123], 

[74] P. Adshead and M. Wyman, Chromo-Natural Inflation: Natural inflation on a steep potential 
with classical non-Abelian gauge fields, Phys. Rev. Lett. 108 (2012) 261302, 

[arXiv: 1202.2366], 

[75] A. Maleknejad, M. M. Sheikh-Jabbari, and J. Soda, Gauge Fields and Inflation, Phys. Rept. 
528 (2013) 161-261, [arXiv: 1212.2921], 

[76] J. F. Dufaux, G. N. Felder, L. Kofman, M. Peloso, and D. Podolsky, Preheating with trilinear 
interactions: Tachyonic resonance, JCAP 0607 (2006) 006, [hep-ph/0602144]. 

[77] A. D. Linde, A New Inflationary Universe Scenario: A Possible Solution of the Horizon, 
Flatness, Homogeneity, Isotropy and Primordial Monopole Problems, Phys.Lett. B108 (1982) 
389-393. 


- 31 


[78] Planck Collaboration Collaboration, P. Ade et al., Planck 2013 results. XXII. Constraints 
on inflation, Astron.Astrophys. 571 (2014) A22, [arXiv: 1303.5082]. 

[79] H. Peiris, R. Easther, and R. Flauger, Constraining Monodromy Inflation , JCAP 1309 (2013) 
018, [arXiv: 1303.2616], 

[80] R. Easther and R. Flauger, Planck Constraints on Monodromy Inflation , JCAP 1402 (2014) 
037, [arXiv: 1308.3736], 

[81] M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P. Hertzberg, Oscillons After 
Inflation, Phys.Rev.Lett. 108 (2012) 241302, [arXiv : 1106. 3335]. 

[82] Planck Collaboration Collaboration, P. Ade et al., Planck 2013 Results. XXIV. Constraints 
on primordial non-Gaussianity, Astron.Astrophys. 571 (2014) A24, [arXiv: 1303. 5084]. 

[83] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, New cosmological constraints on 
primordial black holes, Phys.Rev. D81 (2010) 104019, [arXiv: 0912. 5297]. 

[84] C.-M. Lin and K.-W. Ng, Primordial Black Holes from Passive Density Fluctuations, 

Phys.Lett. B718 (2013) 1181-1185, [arXiv: 1206. 1685], 

[85] L. Alabidi, K. Kohri, M. Sasaki, and Y. Sendouda, Observable Spectra of Induced 
Gravitational Waves from Inflation, JCAP 1209 (2012) 017, [arXiv: 1203.4663]. 

[86] Particle Data Group Collaboration, K. Olive et al., Review of Particle Physics, Chin.Phys. 
C38 (2014) 090001. 

[87] L. Kofman, A. D. Linde, and A. A. Starobinsky, Towards the theory of reheating after 
inflation, Phys.Rev. D56 (1997) 3258-3295, [hep-ph/9704452]. 

[88] C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and 
Engineers. 1978. 

[89] J. Martin and D. J. Schwarz, WKB approximation for inflationary cosmological perturbations, 
Phys.Rev. D67 (2003) 083512, [astro-ph/0210090]. 

[90] G. N. Felder and I. Tkachev, LATTICEEASY: A Program for lattice simulations of scalar 
fields in an expanding universe, Comput.Phys.Commun. 178 (2008) 929-932, 
[hep-ph/0011159], 

[91] A. V. Frolov, DEFROST: A New Code for Simulating Preheating after Inflation, JCAP 0811 
(2008) 009, [arXiv: 0809. 4904], 

[92] R. Easther, H. Finkel, and N. Roth, PSpectRe: A Pseudo-Spectral Code for (P)reheating, 
JCAP 1010 (2010) 025, [arXiv: 1005.1921], 

[93] Z. Huang, The Art of Lattice and Gravity Waves from Preheating, Phys.Rev. D83 (2011) 
123509, [arXiv:1102.0227], 

[94] “http://cosmo.kenyon.edu/gabe.html.” 

[95] M. A. Amin, R. Easther, and H. Finkel, Inflaton Fragmentation and Oscillon Formation in 
Three Dimensions, JCAP 1012 (2010) 001, [arXiv: 1009.2505]. 

[96] M. A. Amin, Inflaton fragmentation: Emergence of pseudo-stable inflaton lumps (oscillons) 
after inflation, arXiv: 1006.3075. 

[97] K. D. Lozanov and M. A. Amin, End of inflation, oscillons, and matter-antimatter 
asymmetry, Phys.Rev. D90 (2014), no. 8 083528, [arXiv: 1408.1811]. 

[98] E. I. Sfakianakis, Analysis of Oscillons in the SU(2) Gauged Higgs Model, arXiv: 1210.7568. 

[99] N. Graham, Numerical Simulation of an Electroweak Oscillon, Phys.Rev. D76 (2007) 085017, 
[arXiv: 0706 .4125]. 


- 32 


[100] Planck Collaboration Collaboration, P. Ade et al., Planck 2015. XX. Constraints on 
inflation , arXiv: 1502.02114. 


- 33 


